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\Q ■ ABSTRACT 

We use 73 h of stereoscopic data taken with the MAGIC telescopes to investigate the very high-energy (VHE) gamma-ray emission of the Crab 

>— ' . pulsar. Our data show a highly significant pulsed signal in the energy range from 50 to 400 GeV in both the main pulse (PI) and the interpulse 

' " (P2) phase regions. We provide the widest spectra to date of the VHE components of both peaks, and these spectra extend to the energy range 

' of satellite-borne observatories. The good resolution and background rejection of the stereoscopic MAGIC system allows us to cross-check the 

, correctness of each spectral point of the pulsar by comparison with the corresponding (strong and well-known) Crab nebula flux. The spectra of 

• i-H ■ both PI and P2 are compatible with power laws with photon indices of 4.0 + 0.8 (PI) and 3.42 ± 0.26 (P2), respectively, and the ratio P1/P2 

■ between the photon counts of the two pulses is 0.54 ± 0.12. The VHE emission can be understood as an additional component produced by the 
J_i inverse Compton scattering of secondary and tertiary e ± pairs on IR-UV photons. 
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1. Introduction 

The Crab pulsar is a young ne utron star tha t is the central rem- 
nant of the supernova SN 1054 dMittonll978l) . It is one of the few 
pulsars that have been detected in almost all energies, ranging 
from radio (e.g.. lLvne et alJ fl993) to VHE gamma rays. In the 
highest-energy regime, it was dete cted up to a few tens of GeV 
by Fermi-LAT ( Abdo et al. 2010a), be tween approximately 25 — 
100 G eV by MAGIC (lAliu et al.ll2008tlSaitd2010tlAleksic et alJ 
l20llb and above 100 GeV by VERITAS lAliu et alj|201 lb . The 
light curves and the spectra obtained by these observations sug- 



gest that gamma-ray pulsars have high-altitude emission zones 
that avoid a super-exponential spectral cutoff, which would be 
caused by magnetic pair production. Consequently, the favored 
models to explain the production of gamma rays to at least a few 
GeV are those in which fan-like beams of high-energy electrons 
scan over a large fraction of the outer magnetosph ere, either 
very close to the lig ht cylinder (outer gap model, Chen getal] 
1986; Romani 199 6J) or all along the last open field lines (slot 
gap model. lAronslll983tlMuslimov & Hardingll20041) . 



1 



Aleksic et al.: Phase-resolved Crab pulsar spectra with MAGIC 



For other rotation-powered gamma-ray pulsars beside the 
Crab pulsar, Fermi-LAT observations have shown that their en- 
ergy spectra exhib it exponential cutoffs at around a few GeV 
( Abd o et aLll20 10b). This mild cutoff has been widely accepted 
as a result of the curvature process by e ± migrating along curved 
paths. In this scenario, the cutoff energy corresponds to the high- 
est characteristic curvature-radiation energy of the p articles ac- 
celerated in the magnetosphere (e.g.. lRomanilll996l) . However, 
the spe ctrum of the Crab pu l sar strongly disfa vors an exponential 
cutoff dAleksic et all 1201 it lAliu et al.11201 lh . making this pul- 
sar a counterexample of the general property. Thus, to develop 
pulsar emission theories beyond the widely accepted curvature- 
radiation models, it is essential to examine the detailed phase- 
resolved spectrum of this youngest pulsar in the HE to VHE 
regimes. 



2. Data set and analysis techniques 

The t wo MAGIC telescopes dAleksic et al.1 1201 2t IZanin et all 
1201 lh situated on the island of La Palma (28.8° N, 17.8° W, 
2220 m a. s.l.), use the imaging atmospheric Cherenkov tech- 
nique to detect gamma rays above a few tens of GeVQ. Since 
summer 2009, when the system started operating in stereoscopic 
mode, its background suppression was substantially improved, 
and a sensitivityQ of 0.8 % Crab ne bula units above 250 GeV has 
been achieved ( Aleksic et al .1120121) . 

In the analysis presented here, we used 73 h of good qual- 
ity stereoscopic data from the winter seasons in 2009/2010 and 
2010/2011. Of these data, 43 h w ere taken in the wobble obser- 
vation mode dFomin etalj [l994). and another 30 h in on-source 
observation mode. The data were taken at zenith angles below 
35° to ensure a low threshold. 

For the data anal ysis, we used the standard MAGIC analy - 
sis package MARS dMoraleio et all 120091; lAleksic et alj|2012l) . 
applying the so-called sum cleaninglLombardi et alJ 201 1 ) to 
achieve the lowest possible threshold. For the gamma/hadron 
separation and gamma d irection estimation we apply the random 
forest (RF) technique (Albert et aD l2008h . Because our back- 
ground is dominated by Crab nebula gamma rays instead of 
hadrons already above ~ 120 GeV, we opted for loose and con- 
servative selection cuts. The phase of each event with respect to 
the main radio pulse w as calculated using the TEMP02 pack- 
age (Hobbs et al. 2006) and the monthly e phemerides pub licly 
provided by the Jodrell Bank ObservatorjQ dLvne et al.lll993l) . 

For th e spectra, w e applied the unfolding algorithms de- 
scribed in Albert et al.l d2007l) . This procedure corrects the mi- 
grations and the energy biases expected in the threshold regime. 
During unfolding iterations, the simulated events are reweighted 
each time with the appropriate spectrum derived in the previous 
iteration. 
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3. Results 

3.1. Folded light curves 

We obtained three folded light curves using all data with esti- 
mated energies between 46 - 416 GeV and for two sub-ranges 



1 The nominal threshold in standard trigger mode, defined as the peak 
of the simulated energy distribution for a Crab-nebula-like spectrum 
after all cuts and at low zenith angles, is 75 - 80 GeV. 

2 Defined as the source strength needed to achieve yV cx / JNb] 
50 h effective on-time. 

3 http://www.jb.man.ac.uk/research/pulsar/crab.html 
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Fig. 1. MAGIC folded light curves of the Crab pulsar for our 
total range in estimated energy and for two separate sub-bins. 
The shaded areas are the on-phase regions P1m and P2m (see 
text), the light shaded area is the off-region [0.52 - 0.87]. The 
dashed line is the constant background level calculated from that 
off-region. 



46-138 GeV and 138 -416 GeV (FigfTJi. The median true ener- 
gies of these samples were estimated from simulations to be ap- 
proximately 100 GeV, 80 GeV and 180 GeV, respectively. The 
significance of the pulsation was tested with the 7? w test, the H 
test dde Jager et al.lfT989h . and a simple ^- 2 -test. None of these 
tests makes an a priori assumption concerning the position and 
the shape of the pulsed emission, and they yield significances 
of 8.6 cr, 6.4 cr and 7.7 cr, respectively. The folded light curve 
clearly shows two distinct peaks, the well-known PI and P2. 

We fitted a very fine-binned all-energy folded light curve, 
maximizing a Poissonian likelihood function that includes two 
Gaussians or Lorentzians over a constant background. The fitted 
Gaussian (Lorentzian) positions of PI and P2 are 0.005 + 0.003 
(0.005 + 0.002) and 0.3996 + 0.0014 (0.3993 + 0.0015), respec- 
tively, with corresponding pulse widths (FWHM) of 0.025 + 
0.007 (0.025 + 0.008) and 0.026 + 0.004 (0.023 + 0.004). The 
signal in P2 is strong enough to also be fitted with an asymmet- 
ric Lorentzian, which involves more parameters. The results are 
displayed in Fig. [2] All fits to our data yield very similar likeli- 
hoods, which neither supports nor excludes the presence of the 
tails implied by a Lorentzian function. Furthermore, the asym- 
metric fit does not yield a significant difference in the leading and 
the trailing wings of P2. Hence, we conclude that the conserva- 
tive approach of using a Gaussian parameterization is sufficient 
to describe our peaks. 

Notably, there is a positive excess throughout the region be- 
tween the two peaks. Most prominently, the trailing wing TW1 
= [0.04-0. 14] has an excess corresponding to 3.4 cr in the lower- 
energy bin (46 - 1 38 GeV), which may allow for a significant de- 
tection once more data is collected. A bridge emission between 
the peaks in our lowest-energy light curve is also expe cted if one 
considers that in the Fermi-LAT data presented in Abd o et a"D 
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MAGIC, P1, 46 GeV < E < 416 GeV 
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Fig. 2. Close-up display of the two fitted peaks PI and P2 us- 
ing a finer binning and a smaller range. The blue solid curves 
represent the Gaussian functions that we use to define the signal 
phase intervals. The red dashed curves are the Lorentzian func- 
tions, which allow for wider tails, and the green solid curve is an 
asymmetric Lorentzian function. The latter did not converge for 
PI (see text). 

J2010ah . the bridge emission is evident up to at least 10 GeV, and 
it is denoted as being spectrally harder than the peak emission. 
However, our current significance in the bridge region is too low 
for a spectral analysis and will not be considered in more detail. 

The peaks we found are significantly narrower than those in 
the GeV regime, and along with MAGIC-Mono and VERITAS 
data, a consistent trend from GeV to beyond 100 GeV can be 
established (Fig. [30. Consequently, the excess we found is 
much more conce ntrated than the wide peak ranges defined in 
iFierroetalJ (119981) (P1 E = [0.94-0.04] andP2 E = [0.32-0.43], 
where E stands for EGRET, in contrast to the MAGIC and 
VERITAS definitions below). Because with too-large phase in- 
tervals one integrates an unnecessarily large number of noise 
events, we decided to investigate the signal both in the EGRET 
intervals and in narrower, a posteriori defined phase in tervals, us- 
ing th e Gaussian peak positions +2 <x, as was done in lAliu et al.1 
(EoTlh . We obtained P1 M = [0.983 - 0.026] and P2 M = [0377- 
0.4221 , the excess of which corresponds to 10.4 cr after iLi & Mai 
d!983[ Eq. 17, see also Table H). The low-/high-energy Li &Ma 
significances for P1 M (P2 M ) are 4.4/3.3 (7.9/5.9). A listing of the 
all-energy significances can be found in Tablefl] 

It is important to note that the two phase interval definitions 
are equally valid. The difference between them is mainly that 
the wide intervals lead to a higher noise contribution but are free 
of any possible selection bias, whereas the narrow intervals have 
much lower noise, but are affected by a mino r selection bias. The 
VERITAS results shown in Aliu et al .1 (1201 11) were calculated us- 
ing Pl v = [0.987 - 0.009] and P2 V = [0.375 - 0.421], which is 
still a bit narrower than our definitions. 

4 A correction to the absolute phase values in Abdo et al. 
(2010a) was announced on the Fermi-LAT websites 
(http://fermi.gsfc.nasa.gOv/ssc/data/access/lat/ephems/l and is in- 
corporated in this plot. 



Crab Pulsar, Pulse extension and phase definitions 
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Fig. 3. Compilation of pulse profi le parameters at di fferent en- 
ergies, measured by F ermi-LAT dAbdo et al.l l2010ai light cir- 
cles), MAGIC-Mono dAleksic et al.l 1201 U sq uares), MAGIC- 
Stereo (this work, dark circles) and VERITAS (lAliu et alJl201ll 
diamonds). The solid points are the phases of the half-maxima 
(PHM), while the crosses indicate the corresponding phases of 
the peak. The vertical lines indicate the phase range definitions 
used for the spectra in Fig. [4] 

For the emission ratio between the two peaks, we found 
0.54 ± 0.12 for P1 M /P2 M and 0.46 ± 0.13 for P1 E /P2 E . We also 
looked for the differences in the pulse shape parameters between 
the two energy intervals of Fig. [T] but we found no significant 
changes in the pulse width, the position, or the relative intensity 
(for either phase range definition). This invariance might be re- 
lated to the fact that although our energy range is almost an order 
of magnitude, the mean energies of the two energy bins (80 vs. 
180 GeV) are comparably close to each other; thus, the lever arm 
is small compared to the energy-dependent trend in Fig. [3] 

3.2. Energy spectra 

We calculated the energy spectra for (P1+P2)m, P1m and P2m, 
which are shown as the red squares in Fig. [4] and for comparison, 
we also calculated the spectra for the unbiased EGRET intervals 
(see above), which are shown as the yellow circles. The latter can 
be compared directly to previous studies, including the mono- 
scopic MAGIC observations. Given that the EGRET intervals 
cover 21 % of the whole phase, they cause a higher background 
noise than the MAGIC phase ranges, which cover only 8.8 %. 
The VERITAS phase intervals cover 6.8 % of the whole phase, 
which is three times less than the EGRET definitions. Although 
most of their narrow pulse may indeed be contained in this inter- 
val, one may expect a certain discrepancy in flux related to this 
difference in selection. 

The spectra we obtain for the EGRET inter vals are com- 
patible with the monoscopic measurements from Aleksic et al. 
(120111) . considering that the statistical deviations are at most ~ 
2 cr and many of the systematic errors of the two measurements 
are independent. Our stereoscopic measurements, however, sup- 
port the possibility that the gamma-ray energy of MAGIC- mono 
data m ay have been over-estimated, as already discussed in lSaitol 
(I2010h . 

The EGRET and MAGIC phase definitions do not result in 
significantly different PI and P2 spectra, although the points 
of the latter are systematically somewhat below the former. 
This is self-consistent, because the EGRET intervals enclose the 
MAGIC intervals and shows that the selection bias that affects 
the latter is probably very small. 

To determine the spectr al parameters and x 1 values, we ap- 
plied a forward unfolding dAlbert et al.ll2007l) . which is the most 
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Fig. 4. Compilation of the spectral measurements of MAGIC and 
Fermi-LAT for the two emission peaks PI and P2, separately 
(b, c) and for both peaks together (a). The VERITAS spectrum 
is only available (and shown) for Ply+P2v (light blue squares 
and solid line). For comparison, the Crab nebula measurements 
of MAGIC and Fermi-LAT (excluding the pulsed component) 
are also shown (green solid circles and diamonds, resp.). Points 
of similar color refer to the same phase intervals (dark red for 
MAGIC, yellow for EGRET, blue for VERITAS intervals, and 
green points for the nebula spectrum, see also text). The blueish 
solid line in the upper plot is the model discussed in Sect.|4]and it 
is above the points because it includes the bridge emission. The 
displayed systematic error of the MAGIC-Stereo measurement 
corresponds to a shift of +17 % in E and +19 % in flux. 



robust method to parameterize the data. The spectra could be de- 
scribed by power laws as shown in Table Q] It should be noted 
that the y 1 values that we found are not optimal, especially for 
the spectrum of P1m- However, the significance of this incon- 
sistency (2.6 cr pre-trial) is too low to claim a feature with the 
data we present here, especially if the systematic uncertainties 
are considered. 

The ratio of the normalization constants between PI and 
P2 at 100 GeV is 0.4 + 0.2, which is consistent with the val- 
ues directly derived from the light curves. We cross-checked the 
(P1+P2) M spectrum by comparing the 2009/10 data to 2010/1 1 
data, on- to wobble-mode data, two zenith angle ranges, two 
quality cut levels and four unfolding algorithms, and found that 
the spectrum was stable within the errors. 

To ensure a good understanding of all possible systematic ef- 
fects, we furthermore determined the Crab nebula spectrum from 
the data taken in wobble mode, analyzing the same energy range 
with the same energy binning. The nebula spectrum that we ob- 
tained with our cuts (see Fig.@]i agrees with both the recent Crab 
n ebula analys i s in |Zanin et al.l (1201 ll) and the Fermi-LAT data 
in lAbdo et al.l (l2010al) . which confirms the good performance of 
our spectral analysis down to 46 GeV. Notably, also the Crab 
nebula flux of the lowest-energy point, which is at approximately 
55 GeV, agrees within errors with the function derive d with a 
combined Fermi-LAT/MAGIC fit in lZanin et ail d2011l) . This fit 
function is basically independent of the lowest-energy MAGIC 
point because it is determined by the statistically much more 
precise points at higher and lower energies. From these results, 
we find no indication that the total systematic flux uncertainty is 
beyon d the standard low-energy numbers given in lAleksic et al.l 
(1201 2lF l These systematic uncertainties are 17 % on the energy 
scale and 19 % on the flux normalization, which is displayed in 
the upper panel of Fig. [4] Assuming a photon index of 3.6, the 
total flux uncertainty including a possible energy bias is there- 
fore ~ 44 % at low energies. The uncertainty of the spectral in- 
dex of such a soft spectrum is approximately 0.2. All MAGIC 
spectra shown in Fig. [4] are unfolded; thus, the statistical errors 
are correlated by 20 - 60 %, reflecting our energy resolution and 
bias, which vary fro m 15 —40%, depending on the energy (see 
lAleksic et al.ll2012h . 

Figure|4]also shows t he Fermi-LAT spectr a in the EGRET in- 
tervals as determined in lAleksic et al.l (1201 ll) . They extrapolate 
consistently to the monoscopic and stereoscopic spectra within 
systematic and statistical uncertainties. To estimate a Fermi-LAT 
spectrum for the MAGIC phase range definition, we s ummed up 
the m atching phase-resolved fit functions provided in Abd o et aD 
d20l0ah . Because their flux constants are differential in phase, 
the emission from the partly covered phase intervals could also 
be approximated. We find that our narrower intervals lead to sub- 
stantially lower GeV equivalent flux spectra. 

In general, when comparing our energy spectra to those ex- 
tracted from Fermi-LAT, VERITAS or MAGIC-Mono data, it is 
important to bear in mind that in addition to the different phase 
interval definitions, all of these experiments suffer different and 
energy-dependent systematic uncertainties that may lead to dis- 
crepancies on the order of 10 - 30 % in energy. 



5 This argument may be regarded as a calibration of the pulsar spectra 
on the nebula spectrum, a method that is not applicable on the nebula 
spectrum itself, but holds for any other source. 
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Table 1. Results of the spectral fits. 

Phase Sdct" /looGev 6 Phot, index ^/n.d.f. 1 ' prob. rf 
(P1+P2) M 104 13.0 ± 1.6 3.57 ± 0.27 10374 O04~ 
P1 M 5.5 3.9+ 1.7 4.0 + 0.8 9.3/2 0.01 

P2 M 9.9 8.8 + 1.0 3.42 + 0.26 6.1/5 0.30 
(P1+P2) E 7.7 15.5 + 2.9 3.9 + 0.4 9.5/4 0.05 
P1 E 3.9 6.5 + 2.0 3.3 + 1.0 3.8/2 0.15 

P2 E 8.0 11.2+ 1.9 3.7 + 0.4 7.2/5 0.21 

(fl) Detection significance after Li & Ma ( 1983, Eq. 17) 
<*> Flux at 100 GeV in units of 10" 11 cm- 2 s" 1 TeV" 1 

(c) Number of degrees of freedom taken from the distribution of 
estimated energies, which may deviate from the number of unfolded 
points in Fig. [4] 

(d) The fit probabilities calculated from the x 1 values do not include 
systematic effects. 



4. Discussion and conclusions 

We found a pulsed VHE gamma-ray signal from the Crab pul- 
sar that allows us to present spectra with an unprecedentedly 
broad energy range and phase resolution. For completeness and 
comparison, we pr ovide analy s es for both the previously used 
phase intervals in iFierro etaf] (|T998), and the narrower peaks 
that we find in our folded light curves. The range of our spectra 
is abo ut one order of mag nitude, and, along with t he MAGIC- 
Mono JXieksic et al.ll201 ll) and the Fermi-LAT data (lAbdo et alJ 
1201 Oaf) , comprise the first gamma-ray spectrum of the Crab 
pulsar from 100 MeV to 400 GeV without any gap. On the 
high-energy end, this result agrees with the recently published 
VERITAS spectrum of P1+P2 above 100 GeV, including also 
the positions and the remarkably narrow widths of the two 
pulses. 

To interpret the observed pulsed spectrum in the context of 
the outer-gap (OG) mo del, we follow the same method as de- 
scribed in Sect. 8.2 of Aleks ic et alJ d201 lb . In this framework, 
the VHE compontent of the spectrum is the inverse Compton 
radiation of secondary and tertiary electron pairs on magneto- 
spheric IR-UV photons. To derive the expected gamma-ray flux 
of this scenario, we solve the set of Poisson equations f or the 
non-corotational potential (Eq. [9] in lAleksic et alJ d201 lb ) with 
the Boltzmann equations for the created electrons and positrons 
and the radiative transfer equation of the emitted photons. 

We present our theoretical calculation of the total pulsed 
spectrum as a violet solid curve in the upper plot of Fig. [4] In 
this calculation, the angle between the rotational and the mag- 
netic axes is assumed to be a — 65°, and the observer's viewing 
angle is £ = 106°. 

In lAleksic et alJd201 ll) . the calculations of both E\\ (the elec- 
tric field component projected along the local magnetic field 
line, which accelerates e*) and the resultant primary gamma-ray 
emissions (curvature+IC) were carried out within 0.7 Rue from 
the rotation axis, where Rue is the radius of the light cylinder. In 
our new calculation, to take account the strong primary IC emis- 
sion that becomes important near the light cylinder, we extend 
the calculation region up to 0.9 Rue for E\\ and up to 1 .5 Rue for 
primary gamma-ray emissions, after confirming that the emis- 
sion above 1.5/?lc is negligible. Here, 0.9 Rue is a safe upper 
boundary for the E\\ calculation, because E\\ is anyway dimin- 
ished at 0.9 Rue owing to the curving-up field-line geometry to- 
wards the rotation axis near the light cylinder. 

A remarkable consequence of this extended calculation is an 
increased inward flux of primary gamma rays originating from 
the upper side of the gap, which leads to a higher abundance of 



pair-produced e ± at lower altitudes (< 0.6 Rue). This screens the 
original E\\ and hence reduces the curvature-radiation compo- 
nent in the primary spectrum. This reduction makes our new cal- 
culation more compatible with the Fermi-LAT data at GeV en- 
ergies but does not significantly affect the secondary and tertiary 
components at energies beyond a few tens of GeV. Hence, we 
conclude that our revised model can reproduce t he total pulsed 
spectr um between 1 and 400 GeV well (see also iLvutikov et all 
j2"0T!t) for an analytical argument of this process). 

However, a remaining caveat of our new calculation is that 
it still includes the bridge emission that is not contained in the 
spectra of only PI +P2. Therefore, it is still above the Fermi-LAT 
flux points in Fig. [4] A phase-resolved modeling is ongoing and 
will be presented in the future. In general, it is however difficult 
to compute the spectral shape above 100 GeV with high preci- 
sion in the present OG model for the Crab pulsar. This is because 
the photon-photon cross section, and therefore the gamma-ray 
absorption, depends on the square of the collision angle, which 
is typically a few degrees. Hence a small variation in the geom- 
etry can have a large impact on the flux that escapes the pulsar. 
Thus, our model should not be interpreted as a hard quantata- 
tive prediction; instead, it is meant to show that the hard compo- 
nent we see in the experiment can quantitatively be met within 
the present understanding of the OG model. Similarly, the slight 
modulations of the power law component are not to be inter- 
preted as a significantly predicted feature. 

Other possible Ansatze to explain the VHE emission in- 
clude the production of inverse Compton radiation in the un- 
shocked pulsar wind outside the light cylinder by pulse d photons 
dAharonian et alj2012l;lAharonian & Bogovalovl2 003). a striped 
pulsar w ind (Petri 201 1), or the annular gap model presented in 
iDu et alJ d2012l) . The two crucial spectral features to establish to 
test these models are the expected spectral upward-kink in the 
transition region between the curvature and the hard component, 
and the detection or exclusion of a terminal cutoff at a few hun- 
dred GeV. 

Another topic that we will be able to address with a 2 - 3 
times larger dataset is the energy dependence of the pulse shape 
parameters. The narrowness of the pulses and its evolution with 
energy are a stringent requirement that the theoretical model- 
ing must fulfill because the folded light curve is almost sta- 
ble against systematic uncertainties. Moreover, the indication of 
pulsed emission in the trailing wing of PI may indicate that a 
VHE signal between the two peaks might be within reach for 
low-threshold IACT systems. The MAGIC telescopes, which are 
being upgraded in 201 1/12, can address these topics in the com- 
ing years when more data will improve the statistical precision 
of the measurements. 
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